Laboratorio #1

Autor. Felipe Urrutia Vargas

Parte A

Uso de $\texttt{Python}$

Ejercicio 1

Ecuaciones diferenciales ordinarias. Dado el sistema de ecuaciones diferenciales no lineal \begin{cases} \dot{x} = 2x − xy, \\ \dot{y}= −y + cos(y), \end{cases} con condiciones iniciales $x_0 = 1$ e $y_0 = 5$.

a)

Encuentre la solución general del sistema sin considerar condiciones iniciales en $\texttt{Python}$ (use $\texttt{dsolve}$ de $\texttt{sympy}$)

In [1]:
!pip install sympy

Requirement already satisfied: sympy in c:\users\felip\anaconda3\lib\site-packages (1.10.1)
Requirement already satisfied: mpmath>=0.19 in c:\users\felip\anaconda3\lib\site-packages (from sympy) (1.2.1)
In [2]:
from sympy import *
In [3]:
%%time
# Definimos las variables
t = symbols("t")
x = Function("x")
y = Function("y")
# Definimos las ecuaciones
eq = [
    Eq(x(t).diff(t), 2*x(t) - x(t)*y(t)), 
    Eq(y(t).diff(t), -y(t) + cos(y(t)))
]
# Solucionamos el sistema
sol = dsolve(eq, funcs=[x(t), y(t)], t=t)
sol

Wall time: 4.72 s
Out[3]:
[Eq(x(t), exp(-C1 + exp(-C2 - t) + 2)),
 Eq(y(t), C1 + log(exp(-C1 + exp(-C2 - t) + 2)))]

Solución general del sistema sin considerar condiciones iniciales

\begin{cases} x(t) &= \exp(-C1 + \exp(-C2 - t) + 2))\\ y(t) &= C1 + \log(\exp(-C1 + exp(-C2 - t) + 2)))\\ \end{cases}

b)

Resuelva el sistema de forma numérica en el intervalo de tiempo $[0, 5]$ en $\texttt{Python}$ mediante $\texttt{solve_ivp}$ de $\texttt{scipy}$. Grafique las soluciones.

In [4]:
import numpy as np
from scipy.integrate import solve_ivp

Definimos el sistema

In [5]:
def OED_system(t, state):
    x, y = state
    dx = 2*x - y*x
    dy = -y+np.cos(y)
    return [dx, dy]
In [6]:
%%time
# Condicion inicial
icv = [1, 5]
# Intervalo de tiempo
t_span = (0.0, 5.0)
t = np.arange(t_span[0], t_span[1], 0.1)
# Resolver
sol_ivp = solve_ivp(OED_system, t_span, icv, dense_output=True)
z = sol_ivp.sol(t)

Wall time: 0 ns
In [7]:
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
In [8]:
fig = go.Figure(data=[
    go.Scatter3d(
        x=z[0, :], 
        y=z[1, :], 
        z=t,
        mode='markers',
        marker=dict(
        size=8,
        color=t,                # set color to an array/list of desired values
        colorscale='Viridis',   # choose a colorscale
        opacity=1
    ))
])
fig.update_layout(
    width = 800,
    height = 500,
    title = "Solución del sistema (z=Tiempo)"
)
fig.show(renderer="notebook")

c)

Utilice el applet $\texttt{pyplane}$ (disponible en https://github.com/TUD-RST/pyplane) para dibujar los diagramas de fase del sistema para distintas condiciones iniciales. En particular muestre el mismo punto inicial utilizado antes.

In [9]:
import run_pyplane

Estudiaremos cuatro condiciones iniciales diferentes, incluyendo la del enunciado.

  • (1, 5)
  • (-1, 5)
  • (1, 0.8)
  • (5, 3)

Condiciones iniciales $(x_0, y_0) = (1, 5)$

Condiciones iniciales $(x_0, y_0) = (-1, 5)$

Condiciones iniciales $(x_0, y_0) = (1, 0.8)$

Condiciones iniciales $(x_0, y_0) = (5, 3)$

Ejercicio 2

Considere el sistema de control en $\mathbb{R}^2$ $$ \dot{x}(t) = Ax(t) + Bu(t) \quad \text{con} \quad A= \begin{bmatrix} 2 & 2 \\ -5 & 3\end{bmatrix}, \quad B=\begin{bmatrix} 5 \\ 1\end{bmatrix} $$ con condiciones indiciales $x_0= (0,0)'$. Resuelva el sistema de forma numerica en el intervalo de tiempo $[0,10]$ en $\texttt{Python}$ mediante $\texttt{solver-ivp}$ para los siguientes controles de lazo abierto (open loop): $$ u(t) = 0.2; \quad u(t) = e^t; u(t) = e^{-t}cos(t). $$ Grafique las soluciones.

In [10]:
# Definimos el sistema asociado al problema de control optimo
def control_system(t, x, u, A, B):
    dx = A @ x + B * u(t)
    return dx
In [11]:
# Definimos los datos del problema
A = np.array([
    [2, 2], 
    [-5, 3]
])
B = np.array([
    5,
    1
])

# Definimos los controles
def u1(t):
    return 0.2

def u2(t):
    return np.exp(t)

def u3(t):
    return np.exp(-t)*np.cos(t)
# Definimos la grilla de los datos
icv = np.array([0, 0])
t_span = (0.0, 10.0)
t = np.arange(t_span[0], t_span[1], 0.01)

Resolvemos el sistema para cada control

In [12]:
%%time
sol_ivp1 = solve_ivp(control_system, t_span, icv, args=[u1, A, B], dense_output=True)
z1 = sol_ivp1.sol(t)

Wall time: 0 ns
In [13]:
%%time
sol_ivp2 = solve_ivp(control_system, t_span, icv, args=[u2, A, B], dense_output=True)
z2 = sol_ivp2.sol(t)

Wall time: 0 ns
In [14]:
%%time
sol_ivp3 = solve_ivp(control_system, t_span, icv, args=[u3, A, B], dense_output=True)
z3 = sol_ivp3.sol(t)

Wall time: 15.6 ms
In [15]:
def plot_sol(t, z, legend="", i=1):
    fig = go.Figure(data=[
        go.Scatter3d(
            x=z[0, :], 
            y=z[1, :], 
            z=t,
            mode='markers',
            marker=dict(
            size=8,
            color=t,                # set color to an array/list of desired values
            colorscale='Viridis',   # choose a colorscale
            opacity=1
        ))
    ])
    fig.update_layout(
        width = 800,
        height = 500,
        title = f"Solución del sistema {legend} (z=Tiempo)"
    )
    fig.show(renderer="notebook")

Graficamos las soluciones

In [16]:
plot_sol(t, z1, legend="para el control u1", i=1)
In [17]:
plot_sol(t, z2, legend="para el control u2", i=2)
In [18]:
plot_sol(t, z3, legend="para el control u3", i=3)

Ejercicio 3

Optimizacion lineal. Resuelva con $\texttt{Python}$ el siguiente problema de programacion lineal (PL). Utilice para ello $\texttt{linprog}$ de $\texttt{scipy}$ y explore los distintos metodos disponibles. \begin{cases} \min &f(x, y, z) = -8x-y-3z\\ \text{s.a.} & -x+y+z \leq 13 \\ & 3x+5y+5z \leq 10 \\ & 9x-5y\leq 20 \\ & x, y, z \geq 0 \end{cases}

In [19]:
from scipy.optimize import linprog
In [20]:
# Costo
c = np.array([-8, -1, -3])
# Restricciones
## Desigualdad
A = np.array([
    [-1, 1, 1], 
    [3, 5, 5],
    [9, -5, 0]
])
b = np.array([13, 10, 20])
## Igualdad
E = None
e = None
## Signo
bnd = [(0, None) for _ in range(c.shape[0])]
In [21]:
def solve_linprog(A, b, E, e, bnd, method="simplex"):
    opt = linprog(c=c, A_ub=A, b_ub=b,
                      A_eq=E, b_eq=e, bounds=bnd,
                      method=method)
    sol = opt.x
    return sol

Calculamos los tiempos para los distintos metodos disponibles

In [22]:
%%timeit
solve_linprog(A, b, E, e, bnd, method = "simplex")

2.38 ms ± 209 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [23]:
%%timeit
solve_linprog(A, b, E, e, bnd, method = "revised simplex")

1.66 ms ± 41 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [24]:
%%timeit
solve_linprog(A, b, E, e, bnd, method = "interior-point")

3.62 ms ± 209 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [25]:
%%timeit
solve_linprog(A, b, E, e, bnd, method = "highs-ipm")

1.14 ms ± 54.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [26]:
%%timeit
solve_linprog(A, b, E, e, bnd, method = "highs-ds")

692 µs ± 13.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [27]:
%%timeit
solve_linprog(A, b, E, e, bnd, method = "highs")

684 µs ± 20.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Calculamos la solucion para el metodo más rapido

In [28]:
sol = solve_linprog(A, b, E, e, bnd, method = "highs")
sol
Out[28]:
array([2.5, 0.5, 0. ])

Los métodos desde el más lento al más rapido:

  1. interior-point
  2. simplex
  3. revised simplex
  4. highs-ipm
  5. highs-ds / highs

Solución:

  • $x=2.5$,
  • $y=0.5$,
  • $z=0.0$

Ejercicio 4

Optimizacion no lineal. Un paquete postal (en forma de paralelepípedo rectangular) debe satisfacer que su altura más el perímetro de su base no puede exceder los 53 cm. Se pretende diseñar un paquete tal que cumpla con esta restricción y que además posea un volumen máximo. Escriba un modelo matemático para este problema y resuélvalo utilizando el comando minimize de scipy.

Figura 1.e4. Paralelepipedo P(h, a, b) de altura h y lados de la base (a, b).

Sea $P(h, a, b)$ un paralelepipedo como en la figura, entonces su volumen $V(h, a, b)$ es igual a $abh$. La suma del perimetro de la base más la altura es $2 (a + b) + h$. Con esto, el problema de optimizacion para encontrar $(h, a, b)$ optimos es equivalente a resolver el siguiente problema no lineal:

\begin{cases} \min & - V(h, a, b) = - abh \\ \text{s.a.} & h + 2(a + b) \leq 53\\ & h, a, b \geq 0 \end{cases}
In [29]:
from scipy.optimize import minimize
from scipy.optimize import Bounds
from scipy.optimize import LinearConstraint
In [30]:
# Funcion objetivo
fun = lambda x: -np.prod(x)
# Punto de partida
x0 = np.ones(3)
# Restricciones del problema
## Lineales
linear_constraint = LinearConstraint([[1, 2, 2]], [0], [53])
## Signo
bounds = Bounds(np.zeros(3), np.inf*np.ones(3))
In [31]:
%%time
res = minimize(fun, x0, method='SLSQP', bounds=bounds, constraints=[linear_constraint])
print(res.message)
sol = res.x
sol

Optimization terminated successfully
Wall time: 6.01 ms
Out[31]:
array([17.66664669,  8.83333081,  8.83334585])

Solución:

  • $h=17.66664669$,
  • $a=8.83333081$,
  • $b=8.83334585$

Ejercicio 5

Encuentre los puntos de intersección $(x, y)$ de las siguientes cónicas:

  • $2x^2 + y = 1$,
  • $(x − 1/2)^2 − 2(y − 1/4)^2 = 1$.

Para esto, utilice el comando $\texttt{fsolve}$ de $\texttt{scipy}$ y grafique.

In [32]:
from scipy.optimize import fsolve
In [33]:
# Sistema de funciones
def func(x):
    return [2*(x[0]**2) + x[1]-1,
            (x[0]-0.5)**2-2*((x[1]-0.25)**2)-1
           ]

# Root 1
root = fsolve(func, [0, 0])
print("Root 1:", root, func(root))

# Root 2
root = fsolve(func, [-1, 0])
print("Root 2:", root, func(root))

Root 1: [-0.53273178  0.4323937 ] [3.801403636316536e-13, 1.1888268147686176e-12]
Root 2: [-0.82653389 -0.36631654] [-1.962319196024964e-12, 2.553379729874905e-11]

Raices:

  1. $(x, y) = (-0.53273178, 0.4323937)$
  2. $(x, y) = (-0.82653389, -0.36631654)$
In [34]:
from sympy import symbols, plot_implicit, Eq

Graficamos las funciones

In [35]:
x, y = symbols('x y')
plot1 = plot_implicit(Eq(1, 2*(x**2) + y), line_color='steelblue', show=False, title="Grafico de ambas conicas (azul: parabola, rojo: hiperbola)")
plot2 = plot_implicit(Eq(1, (x-0.5)**2-2*((y-0.25)**2)), line_color='crimson', show=False)
plot1.extend(plot2)
plot1.show()

Posee dos puntos de intersección, que coinciden con las dos raices encontradas.

Parte B

Control Optimo: Uso de $\texttt{BOCOP}$

$\texttt{BOCOP}$ es un programa open-source diseñado para resolver problemas de control óptimo tipo Mayer a tiempo final fijo o libre y con restricciones de control y estado. Para instalar el solver de control óptimo BOCOP ver el archivo “bocop.pdf”. Encontrará un manual más detallado en www.bocop.org. BOCOP puede resolver problemas de la forma siguiente

$$ (M) \begin{cases} \min_{u(\cdot)} & J(t_0, y(t_0), t_f, y(t_f)) &(\textit{Criterio})\\ & \dot{y}(t) = f(y(t), u(t)), \forall t \in [t_0, t_f] &(\textit{Dinamica})\\ & \Phi_l\leq \Phi(t_0, y(t_0), t_f, y(t_f)) \leq \Phi_u &(\textit{Condicion de borde})\\ & y_l \leq y(t) \leq y_u, u_l \leq u(t) \leq u_u \forall t \in [t_0, t_f] &(\textit{Cotas})\\ & g_l \leq g(y(t), u(t)) \leq g_u, \forall t \in [t_0, t_f] &(\textit{Restricciones mixtas})\\ \end{cases} $$

donde $y(\cdot) \in \mathbb{R}^n$ es el estado del sistema y $u(\cdot) \in \mathbb{R}^n$ el control. En particular, note que la función objetivo solo depende del tiempo inicial y el final y de los estados iniciales y finales.

Ejercicio 6

Consideremos en este ejercicio el problema de Jackson, que modela las reacciones quimicas de tres componentes $A$, $B$, $C$, $$ A \leftrightarrow^1 B \rightarrow^2 C $$ Las variables de estado son $a$, $b$, $c$, representado las fracciones molares de $A$, $B$, $C$, y $k_1$, $k_2$, $k_3$ son las constantes de velocidad de las reacciones quimicas. Aqui el control $u(t) \in [0,1]$ es la fraccion de catalizador que establece el equilibrio entre las reacciones 1 y 2. El objetivo es el de maximizar la produccion del componente $C$.

El problema de control optimo es el siguiente:

$$ (G) \begin{cases} \max_{u(\cdot)} & c(t_f) \\ & \dot{a}(t) = -u(t)(k_1 a(t) - k_2 b(t)), a(0) = 1,\\ & \dot{b}(t) = u(t)(k_1a(t)-k_2b(t))-(1-u(t))k_3b(t), b(0) = 0,\\ & \dot{c}(t) = (1-u(t))k_3b(t), c(0)=0,\\ &u(t) \in [0,1], \forall t \in [0,t_f]\\ \end{cases} $$

este problema ya esta programado y es uno de los ejemplos que viene con $\texttt{BOCOP}$, el cual se puede encontrar en la carpeta "examples/jackson_basic".

1.

Escriba este problema en la forma $(M)$ e indentifique las funciones $J$, $f$, $\Phi$ y $g$ y las cotas $\Phi_l$, $\Phi_u$, $g_l$, $g_u$, $y_l$, $y_u$, $u_l$ y $u_u$.

Basta considerar las siguientes funciones

  • $J(t_0, y(t_0), t_f, y(t_f)) = -c(t_f)$. Maximizar $c$ es equivalente a menos minimizar $-c$.
  • $f(y(t), u(t)) = \begin{bmatrix} -u(t)(k_1 a(t) - k_2 b(t)) \\ u(t)(k_1a(t)-k_2b(t))-(1-u(t))k_3b(t) \\ (1-u(t))k_3b(t) \end{bmatrix}$
  • $\Phi(t_0, y(t_0), t_f, y(t_f)) = \begin{bmatrix} a(0) \\ b(0) \\ c(0) \end{bmatrix}$
  • $g$ libre

donde $t_0 = 0$.

Las siguientes cotas:

  • $\Phi_l = \Phi_u = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$
  • $g_l, g_u$ libres
  • $y_l = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}, y_u = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}$ fracciones molares
  • $u_l = 0, u_u = 1$.

2.

Haga las modificaciones necesarias para resolver el problema de obtener fracciones molares iguales a: $a(t_f)=0.7$, $0.9$, con $t_f \in \{5, 7, 10\}$. Grafique la solución y el control óptimo.

Basta modificar $\Phi$:

  • $\Phi(t_0, y(t_0), t_f, y(t_f)) = \begin{bmatrix} a(0) \\ b(0) \\ c(0) \\ a(t_f)\end{bmatrix}$

y sus cotas respectivas:

  • $\Phi_l = \Phi_u = \begin{bmatrix} 1 \\ 0 \\ 0 \\ a_f \end{bmatrix}$

con $a_f = 0.7, 0.9$ para $t_f \in \{5, 7, 10\}$.

Graficamos las soluciones y el control

In [36]:
# https://github.com/johnny-godoy/bocop-reader
from bocop_solution import BOCOPSolution
In [37]:
import os
import sys

import matplotlib.pyplot as plt
import numpy as np
In [38]:
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)


plt.style.use('ggplot')
plt.rc('axes', titlesize=14)
plt.rc('legend', fontsize=14)
plt.rc('xtick', labelsize=12)
plt.rc('ytick', labelsize=12)
plt.rcParams.update({'font.size': 16})
plt.rcParams['axes.titlesize'] = 16
plt.rcParams["figure.figsize"] = (10, 6)
plt.rcParams.update({'lines.markeredgewidth': 1})
plt.rcParams.update({'errorbar.capsize': 2})

Cargamos las soluciones y controles

In [39]:
file_sols_p2 = {
    "tf: 5, a(tf): 0.7": "D:/github/MA4703/lab1/ex6/a_tf_0_7_tf_5",
    "tf: 5, a(tf): 0.9": "D:/github/MA4703/lab1/ex6/a_tf_0_9_tf_5",
    "tf: 7, a(tf): 0.7": "D:/github/MA4703/lab1/ex6/a_tf_0_7_tf_7",
    "tf: 7, a(tf): 0.9": "D:/github/MA4703/lab1/ex6/a_tf_0_9_tf_7",
    "tf: 10, a(tf): 0.7": "D:/github/MA4703/lab1/ex6/a_tf_0_7_tf_10",
    "tf: 10, a(tf): 0.9": "D:/github/MA4703/lab1/ex6/a_tf_0_9_tf_10",
}
pc = False
sols_p2 = {
    k: BOCOPSolution(v) if pc else BOCOPSolution(v.replace("D:/github/MA4703/lab1/", ""))
    for k, v in file_sols_p2.items()
}
In [40]:
def plot_sol_by_param(sols, params):
    N = len(params)
    fig, ax = plt.subplots(N, 4, figsize=(18, 5*N))
    for i, param in enumerate(params):
        sol = sols[param]
        controls = sol.controls
        states = sol.states
        states["a"].plot(ax[i, 0])
        states["b"].plot(ax[i, 1])
        states["c"].plot(ax[i, 2])
        controls["u"].plot(ax[i, 3]);
        ax[i, 0].set_ylabel(param)
        ax[i, 1].set_ylabel("") 
        ax[i, 2].set_ylabel("")
        ax[i, 3].set_ylabel("")
        if i!=N-1:
            if i==0:
                ax[i, 0].set_title("variable: a")
                ax[i, 1].set_title("variable: b")
                ax[i, 2].set_title("variable: c")
                ax[i, 3].set_title("control: u")

            ax[i, 0].set_xlabel("")
            ax[i, 1].set_xlabel("");
            ax[i, 2].set_xlabel("");
            ax[i, 3].set_xlabel("");

Soluciones con $a(t_f) = 0.7$ y $a(t_f) = 0.9$ para tiempo tf = 5

In [41]:
plot_sol_by_param(sols=sols_p2, params=["tf: 5, a(tf): 0.7", "tf: 5, a(tf): 0.9"])

Soluciones con $a(t_f) = 0.7$ y $a(t_f) = 0.9$ para tiempo tf = 7

In [42]:
plot_sol_by_param(sols=sols_p2, params=["tf: 7, a(tf): 0.7", "tf: 7, a(tf): 0.9"])

Soluciones con $a(t_f) = 0.7$ y $a(t_f) = 0.9$ para tiempo tf = 10

In [43]:
plot_sol_by_param(sols=sols_p2, params=["tf: 10, a(tf): 0.7", "tf: 10, a(tf): 0.9"])

3.

Ahora estamos interesados en maximizar la producción de $C$ pero con cantidades molares, al final del tiempo para $A$ y $B$, por encima de $0.8$ y $0.05$ respectivamente, para cada uno de los $t_f$ del item anterior. Haga las modificaciones necesarias para resolver este nuevo problema. Grafique la solución y el control óptimo.

Cargamos las soluciones y controles

In [44]:
file_sols_p3 = {
    "tf: 5, a(tf)>0.8, b(tf)>0.05": "D:/github/MA4703/lab1/ex6/a_0_8_b_0_05_tf_5",
    "tf: 7, a(tf)>0.8, b(tf)>0.05": "D:/github/MA4703/lab1/ex6/a_0_8_b_0_05_tf_7",
    "tf: 10, a(tf)>0.8, b(tf)>0.05": "D:/github/MA4703/lab1/ex6/a_0_8_b_0_05_tf_10",
}
pc = False
sols_p3 = {
    k: BOCOPSolution(v) if pc else BOCOPSolution(v.replace("D:/github/MA4703/lab1/", ""))
    for k, v in file_sols_p3.items()
}

Soluciones con $a(tf)>0.8$, $b(tf)>0.05$ para tiempos tf = 5, 7, 10

In [45]:
plot_sol_by_param(sols=sols_p3, params=sols_p3.keys())

Ejercicio 7

Considere el siguiente problema de control óptimo:

$$ \min_{u(\cdot)} \int_{0}^{t_f} u(t)^2dt; \quad y^{(3)}(t) = u(t); \quad y(t)\geq 0, $$

con $t_f = 10$, $y(0)=1$, $\dot{y}(0)=-2$, $\ddot{y}(0) = 0$, $\dot{y}(t_f) = 0$ y $\ddot{y}(t_f) = 0$. Note que el criterio depende del control entonces se debe modificar el problema para resolver con $\texttt{BOCOP}$.

Se recomienda revisar una variantes de este problema que ya esta programado y es uno de los ejemplos que viene que con $\texttt{BOCOP}$, el cual se puede encontrar en la carpeta "examples/robbins",

1.

Primeramente identifique el sistema que modela este problema y luego introduzca una nueva variable de estado $z$ tal que $z(t_f ) = \int^{t_f}_0 u(t)^2 dt$. ¿Cuál es la dinámica y la condición inicial de $z$?

El sistema que modela este problema es el siguiente:

$$ (B) \begin{cases} \min_{u(\cdot)} & \int_{0}^{t_f} u(t)^2dt \\ & \dddot{y}(t) = u(t)\\ & y(t) \geq 0\\ & y(0)=1\\ & \dot{y}(0)=-2\\ & \ddot{y}(0) = 0\\ & \dot{y}(t_f) = 0\\ & \ddot{y}(t_f) = 0\\ & u(t) \in \mathbb{R}, \forall t \in [0,t_f]\\ & t_f = 10\\ \end{cases} $$

Ahora bien consideremos $z(s) = \int^{s}_0 u(t)^2 dt$. Con la siguiente dinamica: $$ \dot{z}(s) = u(s)^2, \forall s \in [0,t_f], \quad z(0) = 0 $$ En particular la dinamica de $z(\cdot)$ satisface que $z(t_f) = \int^{t_f}_0 u(t)^2 dt$.

2.

Usando esta variable adicional, escribe un problema de Mayer (M) que es equivalente al problema (B). Identifique las funciones $J$, $f$, $\Phi$ y $g$.

El problema de Mayer ($M$) equivalente a ($B$) esta dado como sigue: $$ (M) \begin{cases} \min_{u(\cdot)} & z(t_f)\\ & \dot{z}(t) = u(t)^2\\ & \dot{c}(t) = u(t)\\ & \dot{b}(t) = c(t)\\ & \dot{a}(t) = b(t)\\ & a(t) \geq 0\\ & z(0) = 0\\ & a(0)=1\\ & b(0)=-2\\ & c(0) = 0\\ & b(t_f) = 0\\ & c(t_f) = 0\\ & u(t) \in \mathbb{R}, \forall t \in [0,t_f]\\ & t_f = 10\\ \end{cases} $$ En este caso las funciones vienen dadas por: $y_M(\cdot) = (z, c, b, a)(\cdot)$ y $y_B(\cdot) = a(\cdot), \dot{y_B}(\cdot) = b(\cdot), \ddot{y_B}(\cdot) = c(\cdot)$.

  • $J(t_0, y(t_0), t_f, y(t_f)) = z(t_f)$
  • $f(y(t), u(t)) = \begin{bmatrix} u(t)^2 \\ u(t) \\ c(t) \\ b(t) \end{bmatrix}$
  • $\Phi(t_0, y(t_0), t_f, y(t_f)) = \begin{bmatrix} z(0) \\ c(0) \\c(t_f)\\ b(0) \\ b(t_f) \\ a(0) \end{bmatrix}$
  • $g$ libre

donde $t_0 = 0$.

Las siguientes cotas:

  • $\Phi_l = \Phi_u = \begin{bmatrix} 0 \\ 0 \\0 \\-2\\0\\1\end{bmatrix}$
  • $g_l, g_u$ libres
  • $y_l = \begin{bmatrix} -\infty \\ -\infty \\ -\infty \\ 0 \end{bmatrix}, y_u = \begin{bmatrix} +\infty \\ +\infty \\ +\infty \\ +\infty \end{bmatrix}$. Para $a(t) \geq 0$.
  • $u_l = -\infty, u_u = +\infty$.

3.

Resuelva este problema considerando $y(t) \leq 1$, $5$ en cada instante de tiempo y $\ddot{y}(t_f) = 3$, $5$, $7$. Grafique los resultados encontrados

Hay que cambiar:

  • $y_u = \begin{bmatrix} +\infty \\ +\infty \\ +\infty \\ 1.5 \end{bmatrix}$. Para $a(t) \geq 0$
  • $\Phi_l = \Phi_u = \begin{bmatrix} 0 \\ 0 \\c_f \\-2\\0\\1\end{bmatrix}$

donde $c_f = 3, 5, 7$.

In [46]:
def plot_sol_by_param(sols, params):
    N = len(params)
    fig, ax = plt.subplots(N, 5, figsize=(18+2, 5*N))
    for i, param in enumerate(params):
        sol = sols[param]
        controls = sol.controls
        states = sol.states
        states["crit"].plot(ax[i, 0])
        states["ddy"].plot(ax[i, 1])
        states["dy"].plot(ax[i, 2])
        states["y"].plot(ax[i, 3])
        controls["u"].plot(ax[i, 4]);
        ax[i, 0].set_ylabel(param)
        ax[i, 1].set_ylabel("") 
        ax[i, 2].set_ylabel("")
        ax[i, 3].set_ylabel("")
        ax[i, 4].set_ylabel("")
        if i!=N-1:
            if i==0:
                ax[i, 0].set_title("variable: crit")
                ax[i, 1].set_title("variable: ddy")
                ax[i, 2].set_title("variable: dy")
                ax[i, 3].set_title("variable: y")
                ax[i, 4].set_title("control: u")

            ax[i, 0].set_xlabel("")
            ax[i, 1].set_xlabel("");
            ax[i, 2].set_xlabel("");
            ax[i, 3].set_xlabel("");
            ax[i, 4].set_xlabel("");

Cargamos las soluciones y controles

In [47]:
file_sols_ex_7_p3 = {
    "y<=1, ddy(tf): 3": "D:/github/MA4703/lab1/ex7/y_1_ddy_3",
    "y<=1, ddy(tf): 5": "D:/github/MA4703/lab1/ex7/y_1_ddy_5",
    "y<=1, ddy(tf): 7": "D:/github/MA4703/lab1/ex7/y_1_ddy_7",
    "y<=5, ddy(tf): 3": "D:/github/MA4703/lab1/ex7/y_5_ddy_3",
    "y<=5, ddy(tf): 5": "D:/github/MA4703/lab1/ex7/y_5_ddy_5",
    "y<=5, ddy(tf): 7": "D:/github/MA4703/lab1/ex7/y_5_ddy_7",
}
pc = False
sols_ex_7_p3 = {
    k: BOCOPSolution(v) if pc else BOCOPSolution(v.replace("D:/github/MA4703/lab1/", ""))
    for k, v in file_sols_ex_7_p3.items()
}

Soluciones con $y(t) \leq 1$

In [48]:
plot_sol_by_param(sols=sols_ex_7_p3, params=["y<=1, ddy(tf): 3", "y<=1, ddy(tf): 5", "y<=1, ddy(tf): 7"])

Soluciones con $y(t) \leq 5$

In [49]:
plot_sol_by_param(sols=sols_ex_7_p3, params=["y<=5, ddy(tf): 3", "y<=5, ddy(tf): 5", "y<=5, ddy(tf): 7"])

4.

Considere ahora que se quiere incluir en la función objetivo el termino $α · y$ con $α = 3$. Resuelva el problema para los mismos valores del item anterior.

Ahora consideremos $z(s) = \int^{s}_0 \alpha y(t) + u(t)^2 dt$. Con la siguiente dinamica: $$ \dot{z}(s) = \alpha y(t) + u(s)^2, \forall s \in [0,t_f], \quad z(0) = 0 $$ En particular la dinamica de $z(\cdot)$ satisface que $z(t_f) = \int^{t_f}_0 \alpha y(t) + u(t)^2 dt$.

Cargamos las soluciones y controles

In [50]:
file_sols_ex_7_p4 = {
    "y<=1, ddy(tf): 3": "D:/github/MA4703/lab1/ex7/y_1_ddy_3_a_3",
    "y<=1, ddy(tf): 5": "D:/github/MA4703/lab1/ex7/y_1_ddy_5_a_3",
    "y<=1, ddy(tf): 7": "D:/github/MA4703/lab1/ex7/y_1_ddy_7_a_3",
    "y<=5, ddy(tf): 3": "D:/github/MA4703/lab1/ex7/y_5_ddy_3_a_3",
    "y<=5, ddy(tf): 5": "D:/github/MA4703/lab1/ex7/y_5_ddy_5_a_3",
    "y<=5, ddy(tf): 7": "D:/github/MA4703/lab1/ex7/y_5_ddy_7_a_3",
}
pc = False
sols_ex_7_p4 = {
    k: BOCOPSolution(v) if pc else BOCOPSolution(v.replace("D:/github/MA4703/lab1/", ""))
    for k, v in file_sols_ex_7_p4.items()
}

Soluciones con $y(t) \leq 1$

In [51]:
plot_sol_by_param(sols=sols_ex_7_p4, params=["y<=1, ddy(tf): 3", "y<=1, ddy(tf): 5", "y<=1, ddy(tf): 7"])

Soluciones con $y(t) \leq 5$

In [52]:
plot_sol_by_param(sols=sols_ex_7_p4, params=["y<=5, ddy(tf): 3", "y<=5, ddy(tf): 5", "y<=5, ddy(tf): 7"])